#1.纹理分析方法-----灰度共生矩阵
import cv2
from PIL import Image
import numpy as np
from matplotlib import pyplot as plt
from struct import unpack
from math import floor
from skimage.feature import greycomatrix, greycoprops

#定义 img传进来是一个unit8格式（图像需要为灰度图）
def ml(img):
    nbit=64#定义级数
    window=5#滑窗的长度
    step = [2]#步长
    angle = [0]#角度
    #将图像归一化处理
    img = np.uint8(255.0 * (img - np.min(img))/(np.max(img) - np.min(img)))
    h, w =img.shape
    #调用numpy里面的函数包，进行比较，bins返回索引值，相当于把灰度降级
    bins = np.linspace(0, 256, nbit+1)
    img1 = np.digitize(img, bins) - 1
    # 图像扩充，调用cv2中的函数进行图像扩充，因选择窗口为5，则行列要各增加2行，用到cv2包
    img2 = cv2.copyMakeBorder(img1, floor(window/2), floor(window/2)
                                  , floor(window/2), floor(window/2), cv2.BORDER_REPLICATE) 
    #用来装各个选的的区域
    patch = np.zeros((window, window, h, w), dtype=np.uint8)
    #对行和列进行循环，在每个像素点进行区域的选定，相当于制定一个滑动的窗口
    for i in range(patch.shape[2]):
            for j in range(patch.shape[3]):
                patch[:, :, i, j] = img2[i : i + window, j : j + window]
    #计算灰度共生矩阵            
    glcm = np.zeros((nbit, nbit, len(step), len(angle), h, w), dtype=np.uint8)
    for i in range(patch.shape[2]):
        for j in range(patch.shape[3]):
            glcm[:, :, :, :, i, j]= greycomatrix(patch[:, :, i, j], step, angle, levels=nbit)
    
    #按公式计算相异性
    def glcm_dissimilarity(glcm, nbit=64):
        dissimilarity = np.zeros((glcm.shape[2], glcm.shape[3]), dtype=np.float32)
        for i in range(nbit):
            for j in range(nbit):
                dissimilarity += glcm[i, j] * np.abs(i-j)
    
        return dissimilarity
    #计算每一个灰度矩阵的相异性，得到的特征会赋值给窗口中心的像素点，则会得到一张特征图片
    for i in range(glcm.shape[2]):        
        for j in range(glcm.shape[3]):
            cut = np.zeros((nbit, nbit, h, w), dtype=np.float32)
            cut = glcm[:, :, i, j, :, :]#转变为4维
            energy1 = glcm_dissimilarity(cut, nbit)
#            也可以直接调用计算的库计算相异性
#            energy1=greycoprops(cut,'dissimilarity')
    return energy1

#读取图片类型，返回数组,一行代表一张图片
def read_image(path):
    with open(path, 'rb') as f:
        magic, num, rows, cols = unpack('>4I', f.read(16))
        img = np.fromfile(f, dtype=np.uint8).reshape(num, 784)
    return img

#读取标签类型，返回数组、
def read_label(path):
    with open(path, 'rb') as f:
        magic, num = unpack('>2I', f.read(8))
        lab = np.fromfile(f, dtype=np.uint8)
    return lab
   
#计算每个图像的相异性
def hb(font):
    save_data=[]
    for i in range(len(font)):
        #提取每一张图片
        eve=font[i].reshape(28,28)
        te=ml(eve).reshape(1,784)
        save_data.append(te)
#        save_data2=np.array(save_data)
        #合并全部数据
        ted=np.concatenate(save_data,axis=0)
    return ted
        
#预测
train=read_image('train-images.idx3-ubyte')
test=read_image('t10k-images.idx3-ubyte')
train_label=read_label('train-labels.idx1-ubyte')
test_label=read_label('t10k-labels.idx1-ubyte')

train_x=hb(train)
test_x=hb(test)
#支持向量机进行预测
from sklearn.svm import LinearSVC
from sklearn.metrics import classification_report
## 初始化线性假设的支持向量机分类器LinearSVC。
lsvc = LinearSVC()

lsvc.fit(train_x, train_label)
#保存模型jj
# 保存模型
#joblib.dump(lsvc, "pima.joblib.dat")

#预测
'''
y_predict = lsvc.predict(test_x)
report=classification_report(test_label,y_predict)                                                                   
'''




#画图表示特征图像
#import matplotlib.pyplot as plt
#plt.figure(figsize=(13, 5))
#font = {'family' : 'Times New Roman',
#'weight' : 'normal',
#'size'   : 12,
#}
#plt.subplot(2,5,1)
#plt.tick_params(labelbottom=False, labelleft=False)
#plt.axis('off')
#plt.imshow(train[0,:].reshape(28,28), cmap ='gray')
#plt.title('Original', font)
#
#
#plt.subplot(2,5,2)
#plt.tick_params(labelbottom=False, labelleft=False)
#plt.axis('off')
#plt.imshow(b, cmap ='gray')
#plt.title('dissimilarity', font)
#
#plt.show()
